library(ggplot2)
library(reshape2)
library(dplyr)
library(tidyr)
library(GGally)
library(grid)
"%&%" = function(a,b) paste(a,b,sep="")
source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan-Paper/scripts/Heather/make-figures/multiplot.R')
my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
fig.dir <- '~/GitHub/GenArch/GenArchPaper/Figures/'
DGN-WB joint heritability. Local h2 is estimated with SNPs within 1 Mb of each gene. Global h2 is estimated with SNPs that are eQTLs in the Framingham Heart Study on other chromosomes (FDR < 0.05).
otherfile<-my.dir %&% 'expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.all.Chr1-22_globalOtherChr.2015-03-18.txt'
fdrother<-read.table(otherfile,header=T) ##FHS eQTLs w/fdr<0.05 on non-gene chromosomes used to define global GRM
##Plot FDR based results
a<-ggplot(fdrother,aes(x=loc.jt.h2,y=glo.jt.h2)) + geom_point(cex=0.8) + geom_abline(intercept=1, slope=-1,color='red') + xlab(expression("local h"^2)) + ylab(expression("global h"^2)) + coord_cartesian(ylim=c(0,1),xlim=c(0,1)) + theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"))
##plot joint h2 estimates
local <- fdrother %>% select(loc.jt.h2,loc.jt.se) %>% arrange(loc.jt.h2)
names(local) = c('h2','se')
data <- local %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE TRUE
## 4974 5989
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE TRUE
## 45.4 54.6
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.13
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))
my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05, y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05, y=0.66, hjust=0,gp=gpar(fontsize=14)))
b<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+coord_cartesian(ylim=c(0,1))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)
global <- fdrother %>% select(glo.jt.h2,glo.jt.se) %>% arrange(glo.jt.h2)
names(global) = c('h2','se')
data <- global %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE TRUE
## 10505 458
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE TRUE
## 95.8 4.2
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.076
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))
my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05, y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05, y=0.66, hjust=0,gp=gpar(fontsize=14)))
c<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("global h"^2)) + xlab(expression("genes ordered by global h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
tiff(filename=fig.dir %&% "Fig1.tiff",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "Fig1.png",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf
## 2
DGN-WB joint heritability with known trans-eQTLs. Local h2 is estimated with SNPs within 1 Mb of each gene. Known trans h2 is estimated with SNPs that are trans-eQTLs in the Framingham Heart Study for each gene (FDR < 0.05).
transfile<-my.dir %&% 'expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.all.Chr1-22_transForGene.2015-03-20.txt'
fdrtrans<-read.table(transfile,header=T) ##FHS trans-eQTLs for gene w/fdr<0.05 used to define Known trans GRM
##Plot FDR based results
a<-ggplot(fdrtrans,aes(x=loc.jt.h2,y=trans.jt.h2)) + geom_point(cex=0.8) + geom_abline(intercept=1, slope=-1,color='red') + xlab(expression("local h"^2)) + ylab(expression("known trans h"^2)) + coord_cartesian(ylim=c(0,1),xlim=c(0,1)) + theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"))
##plot joint h2 estimates
local <- fdrtrans %>% select(loc.jt.h2,loc.jt.se) %>% arrange(loc.jt.h2)
names(local) = c('h2','se')
data <- local %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE TRUE
## 462 732
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE TRUE
## 38.7 61.3
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.133
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))
my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05, y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05, y=0.66, hjust=0,gp=gpar(fontsize=14)))
b<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax,color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+coord_cartesian(ylim=c(0,1),xlim=c(-30,1280))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)
global <- fdrtrans %>% select(trans.jt.h2,trans.jt.se) %>% arrange(trans.jt.h2)
names(global) = c('h2','se')
data <- global %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) )
data <- data[complete.cases(data),]
cigt0 <- data$ymin>0
table(cigt0)
## cigt0
## FALSE TRUE
## 1144 50
ptrue<-round(table(cigt0)/sum(table(cigt0)),3)*100
ptrue
## cigt0
## FALSE TRUE
## 95.8 4.2
meanh2<-round(mean(data$h2),3)
meanh2
## [1] 0.021
data <- mutate(data,nonzeroCI=cigt0,position=1:nrow(data))
my_grob = grobTree(textGrob(substitute(paste("mean ", h^2, " = ", m),list(m=meanh2)), x=0.05, y=0.60, hjust=0,gp=gpar(fontsize=14)))
my_grob2 = grobTree(textGrob(substitute(paste("percent TRUE = ", m),list(m=ptrue[2])), x=0.05, y=0.66, hjust=0,gp=gpar(fontsize=14)))
c<-ggplot(data,aes(x=position,y=h2,ymin=ymin, ymax=ymax,color=nonzeroCI) ) + geom_pointrange(col='gray')+geom_point()+ylab(expression("known trans h"^2)) + xlab(expression("genes ordered by known trans h"^2))+coord_cartesian(ylim=c(0,1),xlim=c(-30,1280))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),legend.justification=c(0,1), legend.position=c(0,1))+annotation_custom(my_grob)+annotation_custom(my_grob2)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
tiff(filename=fig.dir %&% "Fig2.tiff",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "Fig2.png",width=360,height=960)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=1)
dev.off()
## pdf
## 2
Polygenic v. sparse by elastic net.
data<-read.table(my.dir %&% 'working_DGN-WB_exp_10-foldCV_1-reps_elasticNet_eachAlphaR2_hapmap2snps_chr22_2015-01-21.txt',header=T,check.names=F)
ngenes<-dim(data)[1]
print("Elastic Net DGN-WB chr22 (" %&% ngenes %&% " genes)")
## [1] "Elastic Net DGN-WB chr22 (341 genes)"
data_long<-melt(data,by=gene)
## Using gene as id variables
a <- ggplot(data_long, aes(x = as.numeric(levels(variable))[variable] , y = value), group=gene) + geom_line(lwd=0.5,show_guide = FALSE,linetype=1) + aes(color = gene) + xlab(expression(paste("elastic net mixing parameter (",alpha, ")"))) + ylab(expression(paste("10-fold cross-validation R"^2))) + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))+ geom_point(show_guide = FALSE)
print(a)
b<-ggplot(data,aes(x=data[,12],y=data[,2]))+geom_point()+ ylab(expression(paste("R"^2," (",alpha," = 0)"))) + xlab(expression(paste("R"^2," (",alpha," = 0.5)"))) + geom_abline(intercept=0, slope=1,color='red') + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))
c<-ggplot(data,aes(x=data[,12],y=data[,22]))+geom_point()+ ylab(expression(paste("R"^2," (",alpha," = 1)"))) + xlab(expression(paste("R"^2," (",alpha," = 0.5)"))) + geom_abline(intercept=0, slope=1,color='red') + theme_gray(base_size = 18) + coord_cartesian(ylim=c(-0.02,1.02),xlim=c(-0.02,1.02))
data_pairs <- cbind(data[,3],data[,12],data[,21],data[,22])
colnames(data_pairs)<-c("alpha_0.05","alpha_0.5","alpha_0.95","alpha_1")
ggpairs(data_pairs)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)
tiff(filename=fig.dir %&% "Fig3.tiff",width=960,height=320)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "Fig3.png",width=960,height=320)
multiplot(a+ggtitle('A')+theme(plot.title=element_text(hjust=0,face="bold")),b+ggtitle('B')+theme(plot.title=element_text(hjust=0,face="bold")),c+ggtitle('C')+theme(plot.title=element_text(hjust=0,face="bold")),cols=3)
dev.off()
## pdf
## 2
GTEx tissue-wide heritability. Local h2 is estimated with SNPs within 1 Mb of each gene.
h2TW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=T)
seTW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_se.txt",header=T)
###make ordered point+CI h2 plots
gTW_h2 <- gather(h2TW,"Tissue.h2","TissueWide.h2",2:11)
gTW_se <- gather(seTW,"Tissue.se","TissueWide.se",2:11)
gTW_h2_se <- cbind(gTW_h2,gTW_se[,3])
colnames(gTW_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTW_h2_se)/length(unique(gTW_h2_se$Tissue))
gTW_h2_se <- gTW_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0)
fig4<-ggplot(gTW_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),strip.text.x = element_text(size = 18),legend.text=element_text(size=14),legend.title=element_text(size=14))
###calc % nonzero for each tissue TISSUE-WIDE
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-gTW_h2_se %>% select(ensid,Tissue,nonzeroCI) %>% spread(Tissue,nonzeroCI)
for(i in 2:11){
tis<-colnames(a)[i]
cat("\n\n---",tis,"---\n")
print(table(a[,i]))
per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
print(per)
pvec <- c(pvec,per[2])
}
##
##
## --- CrossTissue ---
##
## FALSE TRUE
## 14069 2938
##
## FALSE TRUE
## 82.7 17.3
##
##
## --- AdiposeSubcutaneous ---
##
## FALSE TRUE
## 15880 1127
##
## FALSE TRUE
## 93.40 6.63
##
##
## --- ArteryTibial ---
##
## FALSE TRUE
## 15836 1171
##
## FALSE TRUE
## 93.10 6.89
##
##
## --- HeartLeftVentricle ---
##
## FALSE TRUE
## 16258 749
##
## FALSE TRUE
## 95.6 4.4
##
##
## --- Lung ---
##
## FALSE TRUE
## 16078 929
##
## FALSE TRUE
## 94.50 5.46
##
##
## --- MuscleSkeletal ---
##
## FALSE TRUE
## 15944 1063
##
## FALSE TRUE
## 93.70 6.25
##
##
## --- NerveTibial ---
##
## FALSE TRUE
## 15541 1466
##
## FALSE TRUE
## 91.40 8.62
##
##
## --- SkinSunExposedLowerleg ---
##
## FALSE TRUE
## 15809 1198
##
## FALSE TRUE
## 93.00 7.04
##
##
## --- Thyroid ---
##
## FALSE TRUE
## 15680 1327
##
## FALSE TRUE
## 92.2 7.8
##
##
## --- WholeBlood ---
##
## FALSE TRUE
## 16062 945
##
## FALSE TRUE
## 94.40 5.56
###calc mean h2 for each tissue
a<-gTW_h2_se %>% select(ensid,Tissue,h2) %>% spread(Tissue,h2)
for(i in 2:11){
tis<-colnames(a)[i]
meanh2 <- signif(mean(a[,i]),3)
cat(tis,"mean h2:",meanh2,"\n")
h2vec <- c(h2vec,meanh2)
}
## CrossTissue mean h2: 0.057
## AdiposeSubcutaneous mean h2: 0.0338
## ArteryTibial mean h2: 0.0358
## HeartLeftVentricle mean h2: 0.0383
## Lung mean h2: 0.0325
## MuscleSkeletal mean h2: 0.0279
## NerveTibial mean h2: 0.044
## SkinSunExposedLowerleg mean h2: 0.0351
## Thyroid mean h2: 0.0392
## WholeBlood mean h2: 0.0265
ann_text <- data.frame( h2 = rep(0.78,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
fig4.2<-fig4+geom_text(data=ann_text,aes(label=paste("mean_h^2 ==",mean_h2,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)
ann_text <- data.frame( h2 = rep(0.9,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
fig4<-fig4.2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)
fig4
tiff(filename=fig.dir %&% "Fig4.tiff",width=720,height=960)
fig4
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "Fig4.png",width=720,height=960)
fig4
dev.off()
## pdf
## 2
GTEx cross-tissue and tissue-wide h2 (A) and SE (B).
h2TW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_h2.txt",header=T)
seTW<-read.table(my.dir %&% "GTEx_Tissue-Wide_local_se.txt",header=T)
gh2_TW<-gather(h2TW,"CrossTissue","Tissue",3:11)
colnames(gh2_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5a<-ggplot(gh2_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Wide h'^2)) + ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) + theme(plot.title = element_text(hjust = 0,face="bold"))
gse_TW<-gather(seTW,"CrossTissue","Tissue",3:11)
colnames(gse_TW) <- c('ensid','CrossTissue','TissueName','Tissue')
fig5b<-ggplot(gse_TW,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab('Cross-Tissue SE') + xlab('Tissue-Wide SE') + ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+ theme(plot.title = element_text(hjust = 0,face="bold"))
multiplot(fig5a,fig5b)
tiff(filename=fig.dir %&% "Fig5.tiff",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "Fig5.png",width=480,height=960)
multiplot(fig5a,fig5b)
dev.off()
## pdf
## 2
chr22<-read.table(my.dir %&% "gtex-annot/gencode.v18.genes.patched_contigs.summary.protein.chr22")
gtexEN<-read.table(my.dir %&% "GTEx_tissue-wide_elasticNet_for_ggplot2_2015-05-15.txt",header=T)
chr22EN <- filter(gtexEN,ensid %in% chr22$V5)
ngenes <- length(unique(chr22EN$ensid))
ggplot(chr22EN, aes(x = alpha , y = R2), group=ensid) + facet_wrap(~tissue,ncol=3) + geom_point(show_guide = FALSE) + geom_line(lwd = 0.5, show_guide = FALSE) + aes(color = ensid) + xlab("alpha") + ylab("R2") + ggtitle("GTEx tissue-wide (" %&% ngenes %&% " chr22 genes)")
topR2<-filter(gtexEN,R2>0.5)
ngenes <- length(unique(topR2$ensid))
ggplot(topR2, aes(x = alpha , y = R2), group=ensid) + facet_wrap(~tissue,ncol=3) + geom_point(show_guide = FALSE) + geom_line(lwd = 0.5, show_guide = FALSE) + aes(color = ensid) + xlab("alpha") + ylab("R2") + ggtitle("GTEx tissue-wide (" %&% ngenes %&% " genes with R2>0.5)")
GTEx tissue-specific heritability. Local h2 is estimated with SNPs within 1 Mb of each gene.
h2TS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=T)
seTS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_se.txt",header=T)
###make ordered point+CI h2 plots
gTS_h2 <- gather(h2TS,"Tissue.h2","TissueSpecific.h2",2:11)
gTS_se <- gather(seTS,"Tissue.se","TissueSpecific.se",2:11)
gTS_h2_se <- cbind(gTS_h2,gTS_se[,3])
colnames(gTS_h2_se) <- c("ensid","Tissue","h2","se")
ngenes<-nrow(gTS_h2_se)/length(unique(gTS_h2_se$Tissue))
gTS_h2_se <- gTS_h2_se %>% mutate(ymin = pmax(0, h2 - 2 * se), ymax = pmin(1, h2 + 2 * se) ) %>% arrange(Tissue,h2) %>% mutate(place=rep(1:ngenes,10),nonzeroCI=ymin>0)
figS1<-ggplot(gTS_h2_se,aes(x=place,y=h2,ymin=ymin, ymax=ymax, color=nonzeroCI) ) + facet_wrap(~Tissue,ncol=2) + geom_pointrange(col='gray')+geom_point()+ theme(legend.justification=c(0,1), legend.position=c(0,1))+ylab(expression("local h"^2)) + xlab(expression("genes ordered by local h"^2))+ theme(axis.text=element_text(size=16),axis.title=element_text(size=18,face="bold"),strip.text.x = element_text(size = 18),legend.text=element_text(size=14),legend.title=element_text(size=14))+ coord_cartesian(ylim=c(0,1))
###calc % nonzero for each tissue TISSUE-WIDE
### ADD to plot legend
pvec<-vector()
h2vec<-vector()
a<-gTS_h2_se %>% select(ensid,Tissue,nonzeroCI) %>% spread(Tissue,nonzeroCI)
for(i in 2:11){
tis<-colnames(a)[i]
cat("\n\n---",tis,"---\n")
print(table(a[,i]))
per <- signif(table(a[,i])/sum(table(a[,i])),3)*100
print(per)
pvec <- c(pvec,per[2])
}
##
##
## --- CrossTissue ---
##
## FALSE TRUE
## 14069 2938
##
## FALSE TRUE
## 82.7 17.3
##
##
## --- AdiposeSubcutaneous ---
##
## FALSE TRUE
## 16800 207
##
## FALSE TRUE
## 98.80 1.22
##
##
## --- ArteryTibial ---
##
## FALSE TRUE
## 16701 306
##
## FALSE TRUE
## 98.2 1.8
##
##
## --- HeartLeftVentricle ---
##
## FALSE TRUE
## 16709 298
##
## FALSE TRUE
## 98.20 1.75
##
##
## --- Lung ---
##
## FALSE TRUE
## 16784 223
##
## FALSE TRUE
## 98.70 1.31
##
##
## --- MuscleSkeletal ---
##
## FALSE TRUE
## 16573 434
##
## FALSE TRUE
## 97.40 2.55
##
##
## --- NerveTibial ---
##
## FALSE TRUE
## 16628 379
##
## FALSE TRUE
## 97.80 2.23
##
##
## --- SkinSunExposedLowerleg ---
##
## FALSE TRUE
## 16558 449
##
## FALSE TRUE
## 97.40 2.64
##
##
## --- Thyroid ---
##
## FALSE TRUE
## 16565 442
##
## FALSE TRUE
## 97.4 2.6
##
##
## --- WholeBlood ---
##
## FALSE TRUE
## 16517 490
##
## FALSE TRUE
## 97.10 2.88
###calc mean h2 for each tissue
a<-gTS_h2_se %>% select(ensid,Tissue,h2) %>% spread(Tissue,h2)
for(i in 2:11){
tis<-colnames(a)[i]
meanh2 <- signif(mean(a[,i]),3)
cat(tis,"mean h2:",meanh2,"\n")
h2vec <- c(h2vec,meanh2)
}
## CrossTissue mean h2: 0.057
## AdiposeSubcutaneous mean h2: 0.0204
## ArteryTibial mean h2: 0.0249
## HeartLeftVentricle mean h2: 0.0323
## Lung mean h2: 0.0223
## MuscleSkeletal mean h2: 0.0214
## NerveTibial mean h2: 0.0281
## SkinSunExposedLowerleg mean h2: 0.0288
## Thyroid mean h2: 0.0265
## WholeBlood mean h2: 0.0261
ann_text <- data.frame( h2 = rep(0.78,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
figS1.2<-figS1+geom_text(data=ann_text,aes(label=paste("mean_h^2 ==",mean_h2,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)
ann_text <- data.frame( h2 = rep(0.9,10), place = rep(5000,10), percent= pvec, mean_h2 = h2vec, Tissue = factor(colnames(a)[2:11]), ymin=rep(0.9,10), ymax=rep(0.9,10), nonzeroCI=rep(NA,10), se=rep(0.9,10))
figS1<-figS1.2+geom_text(data=ann_text,aes(label=paste("percent_TRUE ==",percent,sep="")),color="black",show_guide=F,parse=T,hjust=0,size=5)
figS1
tiff(filename=fig.dir %&% "FigS1.tiff",width=720,height=960)
figS1
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "FigS1.png",width=720,height=960)
figS1
dev.off()
## pdf
## 2
GTEx cross-tissue and tissue-specific h2 (A) and SE (B).
h2TS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_h2.txt",header=T)
seTS<-read.table(my.dir %&% "GTEx_Tissue-Specific_local_se.txt",header=T)
gh2_TS<-gather(h2TS,"CrossTissue","Tissue",3:11)
colnames(gh2_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2a<-ggplot(gh2_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab(expression('Cross-Tissue h'^2)) + xlab(expression('Tissue-Specific h'^2)) + ggtitle("A")+ coord_cartesian(ylim=c(-0.1,1.1),xlim=c(-0.1,1.1)) + theme(plot.title = element_text(hjust = 0,face="bold"))
gse_TS<-gather(seTS,"CrossTissue","Tissue",3:11)
colnames(gse_TS) <- c('ensid','CrossTissue','TissueName','Tissue')
figS2b<-ggplot(gse_TS,aes(x=Tissue,y=CrossTissue)) +facet_wrap(~TissueName,scales="fixed",ncol=3) + geom_point(alpha=1/10) + geom_abline(intercept=0, slope=1,color='red') + ylab('Cross-Tissue SE') + xlab('Tissue-Specific SE') + ggtitle("B") + coord_cartesian(ylim=c(-0.01,0.16),xlim=c(-0.01,0.16))+ theme(plot.title = element_text(hjust = 0,face="bold"))
multiplot(figS2a,figS2b)
tiff(filename=fig.dir %&% "FigS2.tiff",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## pdf
## 2
png(filename=fig.dir %&% "FigS2.png",width=480,height=960)
multiplot(figS2a,figS2b)
dev.off()
## pdf
## 2